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Abstract. The present review will focus on recent development of exact-diagonali- 
^ ' zation (ED) methods that use Lanczos algorithm to transform large sparse matri- 

ces onto the tridiagonal form. We begin with a review of basic principles of the 
Lanczos method for computing ground-state static as well as dynamical properties. 
^ ' Next, generalization to finite-temperatures in the form of well established finite- 

temperature Lanczos method is described. The latter allows for the evaluation of 

a temperatures T > static and dynamic quantities within various correlated models. 
Several extensions and modification of the latter method introduced more recently 
, are analysed. In particular, the low-temperature Lanczos method and the micro- 

pH . canonical Lanczos method, especially applicable within the high-T regime. In order 

Q ' to overcome the problems of exponentially growing Hilbert spaces that prevent ED 

O . calculations on larger lattices, different approaches based on Lanczos diagonalization 

within the reduced basis have been developed. In this context, recently developed 
method based on ED within a limited functional space is reviewed. Finally, we briefly 
discuss the real-time evolution of correlated systems far from equilibrium, which can 
be simulated using the ED and Lanczos-based methods, as well as approaches based 
, on the diagonalization in a reduced basis. 

a^ 

1.1 Introduction 

Models of strongly correlated systems have been one of the most intensively 
studied theoretical subjects in the last two decades, stimulated at first by the 
discovery of compounds superconducting at high-temperatures and ever since 
, by the emergence of various novel materials and phenomena which could be 

^ ' traced back to strongly correlated electrons in these systems. Recently, cold 

atoms in optical lattice offer a different realization of strongly correlated quan- 
tum entities, whereby these systems can be even tuned closer to theoretical 
models. 

One of the most straightforward methods to numerically deal with the lat- 
tice (discrete) models of correlated particles, which arc inherently many-body 
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(MB) quantum systems, is exact diagonalization (ED) of small-size systems. In 
view of the absence of well - controlled analytical methods, ED method has 
been employed intensively to obtain results for static and dynamical prop- 
erties of various models with different aims: a) to search and confirm novel 
phenomena specific for strongly correlated systems, b) to test theoretical ideas 
and analytical results, c) to get reference results for more advanced numerical 
techniques. 

MB quantum lattice models of interacting particles are characterized with 
a dimension of the Hilbert space given by the number of basis states Ngt oc 
that is in turn exponentially increasing with the lattice size iV, where K is 
the number of local quantum states. It is therefore clear that ED methods can 
treat fully only systems with limited TV^t, i.e., both K and N must be quite 
modest. 

Among the ED approaches the full ED within the Hilbert space of the 
model Hamiltonian, yielding all eigenenergies and eigenfunctions, is the sim- 
plest to understand, most transparent and easy to implement. In principle it 
allows the evaluation of any ground state (g.s.) property as well as finite tem- 
perature T > static or dynamic quantity, at the expense of very restricted 
Nst- In spite of that, it represents a very instructive approach but also remains 
essentially the only practical method when all exact levels are needed, e.g., 
for studies of level statistics. 

Lanczos-based ED methods have already long history of applications since 
Cornelius Lanczos [1] proposed the diagonalization of sparse matrices using 
the iterative procedure, allowing for much bigger Hilbert spaces N^t relative 
to full ED. Lanczos diagonalization technique is at present a part of standard 
numerical linear algebra procedures [2l[3] and as such in solid state physics 
mainly used to obtain the g.s. energies and wavefunction with corresponding 
expectation values. The approach has been quite early on extended to calcu- 
lation of the dynamical response functions within the g.s. [4]. The method has 
been in the last 20 years extensively used in connection with models related 
to high- Tc materials, for which we can refer to an earlier overview (S). 

The present review will focus on recent development of ED-based and 
Lanczos-based methods. The basics of the Lanczos method are presented in 
Sec ll.2l and its application for g.s. properties in Sec ll.3l One of already estab- 
lished generalizations is the finite-temperature Lanczos method (FTLM) [6ll7] , 
reviewed in Sec ll.4l which allows for the evaluation of T > static and dy- 
namic properties within simplest models. Several extensions and modifica- 
tions of the latter have been introduced more recently, in particular the low- 
temperature Lanczos method (LTLM) [8j and the microcanonical Lanczos 
method (MCLM) [9], particularly applicable within the high-T regime. 

Since the application of the ED methods there have been attempts and 
proposals for the proper choice of reduced basis which could allow for the 
study of bigger systems. While this is clearly very broad subject with most 
substantial achievements in one-dimensional (ID) systems within the frame- 
work of the density-matrix renormalization-group (DMRG) idea, there are 
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also successful applications in higher D> 1 combined with the Lanczos pro- 
cedure being reviewed in Sec ll.51 Recently, there is also quite an intensive 
activity on studies of real-time evolution of correlated systems, both under 
the equilibrium and the non-equilibrium conditions that can be simulated 
using the ED and Lanczos-based methods, as discussed in Sec ll.71 

1.2 Exact diagonalization and Lanczos method 
1.2.1 Models, geometries and system sizes 

ED-based methods are mostly restricted to simple models with only few local 
quantum states K per lattice site in order to reach reasonable system sizes 
. Consequently, there are only few classes of MB models that so far exhaust 
the majority of ED and Lanczos-method studies, clearly also motivated and 
influenced by the challenging physics and relevance to novel materials and 
related experiments. 

To get some feeling for available sizes reachable within the ED-based ap- 
proaches, it should be reminded that in full ED routines the CPU time scales 
with the number of operations Op oc iV^^, while the memory requirement 
is related to the storage of the whole Hamiltonian matrix and all eigenvec- 
tors, i.e., Mem oc N'^i. This limits at present stage of computer facilities the 
full ED method to Nst < 2.10" MB states. On the other hand, using the 
Lanczos-based iterative methods for the diagonalization of sparse matrices 
(Hamiltonians), CPU and memory requirements scale as Op, Mem oc Nst, 
at least in their basic application, to calculate the g.s. and its wavefunction. 
In present-day applications this allows the consideration of much larger basis 
sets, i.e., N^t < 10^. StiU, lattice sizes N reached using the Lanczos technique 
remain rather modest, compared to some other numerical approaches as the 
DMRG and quantum-Monte-Carlo (QMC) methods, if the full Hilbert basis 
space relevant for the model is used. 

The simplest nontrivial class of MB lattice models are spin models, the 
prototype being the anisotropic Heisenberg model for coupled = 1/2 spins. 



where the sum (ij) runs over pairs of lattice sites with an arbitrary interaction 
Jj"" (being in principle anisotropic) and S" are component of local S* = 1/2 
operator. The model has K = 2 quantum states per lattice site and therefore 
allows for biggest N in the ED-based approaches where Nst oc 2^ basis states. 
To reduce Nst as many symmetries and good quantum numbers as practically 
possible are used to decompose the Hamiltonian into separate blocks. Evident 
choice are sectors with the (z-component of) total spin iSf^^ and the wavevector 
q for systems with periodic boundary conditions, but further also rotational 
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symmetries of particular lattices have been used. In this way system sizes up 
to '-^ 36 (for largest and most interesting sector Sf^f = 0) have been reached 
so far using the Lanczos technique without any basis reduction. 

On the basis of this simple model one can already discuss the feasibility 
of the Lanczos-based methods with respect to other numerical quantum MB 
methods. For the g.s. in ID spin systems more powerful methods allowing 
for much bigger systems are DMRG and related approaches. For unfrustrated 
models in D > 1 the QMC methods are superior for the evaluation of static 
quantities at any T. Still, Lanczos-based methods become competitive or at 
least not superseded for frustrated spin models (where QMC can run into 
minus-sign problem) or for dynamical properties at T > 0. 

Next in complexity and very intensively studied prototype model is the 
t-J model, representing strongly correlated itinerant electrons with an anti- 
ferromagnetic (AFM) interaction between their spins. 



where due to the strong on-site repulsion doubly occupied sites are forbidden 
and one is dealing with projected fermion operators Cis = Cis(l — n-i^-s). The 
model can be considered as a good microscopic model for superconducting 
cuprates which are doped Mott insulators. For a theoretical and experimental 
overview of Mott insulators and metal-insulator transitions see Ref. [10] . and 
has been therefore one of the most studied using the Lanczos method [5]. 
It has K = 3 quantum states per lattice site and besides 3^^^. and q, also 
the number of electrons (or more appropriate the number of holes = 
N — iVe) are the simplest quantum numbers to implement. Since the model 
reveals an interesting physics in D > 1, the effort was in connection with 
high- Tc cuprates mostly on 2D square lattice. Here the alternative numerical 
methods have more drawbacks (e.g., minus sign problem in QMC methods 
due to the itinerant character of fcrmions) so Lanczos-based methods are still 
competitive, in particular for getting information on T > dynamics and 
transport. The largest systems considered with Lanczos method so far are 2D 
square lattice with = 32 sites and Nh = 4 holes [TT| . 

Clearly, one of the most investigated within the MB community is the 
standard single-band Hubbard model, which has K = A states per lattice site. 
Due to the complexity Ngt oc 4^ the application of ED and Lanczos-based 
method is already quite restricted reaching so far N = 20 sites [12] requiring 
already Nst ~ 10^ basis states. The model is also the subject of numerous 
studies using more powerful QMC method and various cluster dynamical- 
mean-field-theory (DMFT) methods for much larger lattices so Lanczos-based 
approaches have here more specific goals. 

Since reachable lattices for above mentioned models are rather small it is 
important to choose properly their geometries. This is not the problem in ID 
models, but becomes already essential for 2D lattices, analysed in connection 
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with novel materials, in particular high-Tc cuprates and related materials. In 
order to keep the periodic boundary conditions for 2D square lattice the choice 
of Pythagorean lattices with TV = + with A^; , Xy [T^ has significantly ex- 
tended available sizes. Some of frequently used are presented in Fig ll.ll Taking 
into account only even N such lattices include N = 8, 10, 16, 18, 20, 26, 32, 36 
sites. While unit cells of such lattices are squares, it has been observed that 
they are not always optimal with respect to the number of next-nearest- 
neighbors and further nearest neighbors. It has been claimed and partly tested 
that better result are obtained with slightly deformed lattices (still with peri- 
odic boundary conditions) which at the same time offer an even bigger choice 
of sizes [14]. 
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Fig. 1.1. Tilted clusters used in 2D square- lattice studies 



1.2.2 Lanczos diagonalization technique 

The Lanczos technique is a general procedure to transform and reduce a sym- 
metric Nst X Nst matrix A to a symmetric Mx M tridigonal matrix Tm- From 
the chosen initial A^sf-dimensional vector Vi one generates an orthogonal basis 
of {vi, • • • vm} vectors which span the Krylow space {vi, Avi, • • • A'^^-'-vi} 

In usual applications for the quantum MB system defined with the Hamil- 
tonian operator H the Lanczos algorithm starts with a normalized vector |(/)o), 
chosen as a random vector in the relevant Hilbert space with Ngt basis states. 
The procedure generates orthogonal Lanczos vectors Lm = {\4'm), m = 0, M} 
spanning the Krylow space {|0o), -ff |(/>o) • ■ ■ , H'^^\(j)o)}. Steps are as follows: H 
is applied to |0o) and the resulting vector is split in components parallel to 
\4>o), and normalized \(t>i) orthogonal to it, respectively, 

H\(j,o)^ao\cl)o)+h\cl)i). (1.3) 

Since H is Hermitian, oq = ((/jqI^I^o) is real, while the phase of |0i) can be 
chosen so that bi is also real. In the next step H is applied to |0i). 



i/|<^i)=6'i|0o)+ai|0i)+62|02). 



(1.4) 
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where |02) is orthogonal to {(p^) and \4>i). It foUows also b[ = {(j)Q\H\(f>i) — bi. 
Proceeding with the iteration one gets in i steps 

H\(f>,) = + a, 10,) + b,+,\(f>,+i), l<i< M. (1.5) 

where in Eq. (jl.5l) by construction there are no terms involving \(j)i-2) etc. By 
stopping the iteration at i = Af and setting ^az+i = 0, the Hamiltonian can 
be represented in the basis of orthogonal Lanczos functions \(f>i) as the tridi- 
agonal matrix Hm with diagonal elements ai,i = O.M, and off-diagonal ones 
bi,i — 1,M. Such a matrix is easily diagonalized using standard numerical 
routines to obtain approximate eigenvalues ej and corresponding orthonor- 
mal eigenvectors |V'j); 

M 

IV'j) =^t'r<l<^.): i = 0,M. (1.6) 

4=0 

It is important to realize that are (in general) not exact eigenfunctions 
of H, but show a remainder. On the other hand, it is evident from the diago- 
nalization of Hm that matrix elements 

{i;,\H\i^,) ^e,S,,, t,j = 0,M, (1.7) 

are diagonal independently of Lm (but provided i,j < Af), although the 
values Ej can be only approximate. 

If in the equation ()1.5|) 6m+i — 0, we have found a {M + l)-dimensional 
eigenspace where Hm is already an exact representation of H. This inevitably 
happens when M — Nst — 1, but for M < Nst — 1 it can only occur if the 
starting vector is orthogonal to some invariant subspace of H which we avoid 
by choosing the input vector |(/)o) as a random one. 

It should be recognized that the Lanczos approach is effective only for 
sparse Hamiltonians, characterized by the connectivity of each basis state with 
Kn <C Nst basis states. All prototype discrete tight-binding models discussed 
in Sec ll.2.11 are indeed of such a type in the local MB basis. Estimating the 
computation requirements, the number of operations Op needed to perform 
M Lanczos iterations scales as Op oc KnMNst- The main restriction is still 
in memory requirements due to large Nst- A straightforward application of 
Eq. (|1.5p would require the fast storage of all |0i), « = 0, M, i.e., also the mem- 
ory capacity Mem cx MNst- However, for the evaluation of the eigenvalues 
alone during the iteration, Eq. (|1.5p . only three are successively required, 
so this leads to Mem cx SiV^t . If the Hamiltonian matrix is not evaluated on 
the fly (simultaneously), then also Mem oc K„Nst for the nonzero Hamilton 
matrix elements is needed. 

The Lanczos diagonalization is in essence an iterative power method which 
is known to converge fast for the extreme lower and upper eigenvalues |21[3]) 
clearly in physical application most relevant is the search for the g.s. energy Eq 
and corresponding wavefunction |lfb)- Typically, > 50 are enough to reach 
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very high accuracy for both. It is evident that for such modest M <C N^t one 
cannot expect any reUable results for eigenstates beyond the few at the bot- 
tom and the top of the spectrum. On the other hand, the Lanczos procedure 
is subject to roundoff errors, introduced by the finite-precision arithmetics 
which usually only becomes severe at larger M > 100 after the convergence 
of extreme eigenvalues and is seen as the loss of the orthogonality of vectors 
\4>i)- It can be remedied by successive reorthogonalization [2j^il5] of new 
states |(/)-), plagued with errors, with respect to previous ones. However this 
procedure requires Op ~ M'^Ngt operations, and can become computationally 
more demanding than Lanczos iterations themselves. This effect also prevents 
one to use the Lanczos method, e.g., to efficiently tridiagonalize full large 
matrices 

1.3 Ground State Properties and Dynamics 

After \^q) is obtained, the g.s. static properties can be evaluated in principle 
for any operator A as 

Ao = {%\A\%). (1.8) 

Clearly, the procedure (|1.8p for large basis is effective only if operator A is in 
the same basis also sparse, as it is in most cases of interest. 

It is, however, the advantage of the Lanczos procedure that also g.s. dy- 
namical functions can be calculated within the same framework [4]. Let us 
consider the dynamical (autocorrelation) response function 

C{uj) = (.Z^olAt^-l -Al'Fo), (1.9) 

for the observable given by the operator A where = cj + ie, e > 0. To 
calculate C{uj) one has to run the second Lanczos procedure with a new initial 
function |</)o), 

l^o) = -^Itf'o), a = J{n\A^A\^^). (1.10) 
a ^ 

Starting with \4>q) one generates another Lanczos subspace L^j — j — 

0,M} with (approximate) eigenvectors and eigenenergies tj. The matrix 
for H in the new basis is again a tridiagonal one with a.j and hj elements, 
respectively. Terminating the Lanczos procedure at given M, one can evaluate 
Ea. (|1.9p as a resolvent of the Hj^ matrix expressed in the continued- fraction 
form 
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terminating with bj^jj^-^ — 0, although other termination functions can also be 
employed and well justified. 

We note that frequency moments of the spectral function 

1 r°° 

/ JlmC{uj)dLo = {1'q\A\H - EQfA\^Q) = 

J-oo 

= a^{4>o\iH-Eoy\4>o), (1-12) 

are exact for given \'Pq) provided I < M, since the operator H'',! < M, is 
exactly reproduced within the Lanczos (or corresponding Krylow) space Lj^j. 

Finally, C{lu) (|l.lip can be presented as a sum of j = 0, M poles at 
UJ — ij — Eq with corresponding weights Wj. As a practical matter we note 
that in analogy to Eg. dl.Gp 

w, = 1(^,1^1*^0)1' = «'l(V^j#o)|' = a'v%, (1.13) 

hence no matrix elements need to be evaluated within this approach. In con- 
trast to the autocorrelation function (11.111) . the procedure allows also the 
treatment of general correlation functions Cab(w), with B ^ A"^ . In this case 
matrix elements (!f'o|S|V'j) have to be evaluated explicitly. It should be also 
mentioned that at least lowest poles of C(w), Ea. (|l.lip . should coincide with 
eigenenergies uj — Ei ~ Eq if |0o) is not orthogonal to |!f'o)- However, using 
M > 50 spurious poles can emerge (if no reorthogonalization is used) which, 
however, carry no weight as also evident from exact moments (|1.12p . 

In this chapter we do not intend to present an overview of applications of 
the full ED and Lanczos- type studies of g.s. static and dynamical properties of 
correlated systems. Such investigations have been numerous even before the 
high-Tc era but intensified strongly with studies of prototype models relevant 
for high-T cuprates 'S^ and other novel materials with correlated electrons. Al- 
though variety of models have been investigated they are still quite restricted 
in the number of local degrees and sizes. 



1.4 Static Properties and Dynamics at T > 

Before describing the Finite temperature Lanczos method (FTLM) we should 
note that the Lanczos basis is very useful and natural basis to evaluate the 
matrix elements of the type 

Wki = {n\H''BH'A\n), (1.14) 

where \n) is an arbitrary normalized vector, and A, B are general operators. 
One can calculate this expression exactly by performing two Lanczos pro- 
cedures with M = max(A;, I) steps. The first one, starting with the vector 
\4>o) = W), produces the Lanczos basis Lm along with approximate eigenstates 
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and ej. The second Lanczos procedure is started with the normahzed vec- 
tor |0o) ^If/io) = ^l"-), Ea. (|1.10p . and generates Lm with corresponding jV'j) 
and tj . We can now define projectors onto hmited subspaces 

M M 

Fm = ^|V'.)(V'.|, PM-^|Vi4)(Vi»|. (1-15) 

1=0 i=0 

Provided that (/, k) < M projectors Pm and Pm span the whole relevant 
basis for the operators H'' and i?', respectively, so that one can rewrite Wki 
in Eg. dTTTil) as 

Wki = {MPmHPmH ...HPmBPmH ...PMHPMA\<j,o). (1.16) 

Since H is diagonal in the basis jV'j) and respectively, one can write 

finally 

M M 
i=0 j=0 

It is important to note that matrix element expression (|1.17p is exact, inde- 
pendent of how (in)accurate representation and respectively, 
are to true system eigenvalues. The only condition remains that number of 
Lanczos steps is sufficient, i.e., M > (Z,fc). 

1.4.1 Finite temperature Lanczos method: Static quantities 

A straightforward calculation of canonical thermodynamic average of an op- 
erator A at T > (in a finite system) requires the knowledge of all eigenstates 
\^n) and corresponding energies Em obtained, e.g., by the full ED of i?, 

(A) =5^e-^^"(^„|A|f„) / ^e-^^", (1.18) 

n—l ' n—1 

where (3 — l/ksT. Such direct evaluation is both CPU time and storage 
demanding for larger systems and is at present accessible only for Ngt ~ 10000. 

In a general orthonormal basis |n) for finite system with N^t basis states 
one can express the canonical expectation value (A) as 

(A) = J2{n\e-'''A\n) / (1.19) 

n—l ' n—l 

The FTLM for T > is based on the evaluation of the expectation value in 
Eq. (|1.19p for each starting \n) using the Lanczos basis. We note that such 
procedure guarantees correct high-T expansion series (for given finite system) 
to high order. Let us perform the high-T expansion of Eg. (II. 191) . 
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nst OO (_n\k 

{A) = z-^Y.Y.-^^^\H'm. 

n=l k=a 
n=l k=0 

Terms in the expansion {n\H^A\n) can be calculated exactly using the Lanczos 
procedure with M > k steps (with |0q) — \n) as the starting function) since 
this is a special case of the expression (|1.14l) . Using relation (|1.17p with I = 
and _B = 1, we get 

M 

{n\H''A\n} = (Cl^l") (^D'- (1-21) 

4=0 

Working in a restricted basis k < M ,vfe can insert the expression (|1.21l) into 
sums (jl.20p . extending them to fc > M. The final result can be expressed as 

n=l 1=0 

N,t M 

^-EE^"'"'' wrx^ri"), (1.22) 

n=l i=0 

and the error of the approximation is 0(/3^^+^). 

Evidently, within a finite system Eq. (|1.22p . expanded as a series in /?, 
reproduces exactly the high-T series to the order M . In addition, in contrast 
to the usual high-T expansion, Eq. (|1.22p remains accurate also for T 0. Let 
us assume for simplicity that the g.s. jS'o) is nondegenerate. For initial states 
\n) not orthogonal to |!f'o), already at modest M ^ 50 the lowest eigenstate 
IV'o) converges to We thus have for /3 — > oo, 

{A) = Y,{n\%){']^o\A\n) / ^(n|.Z'o)('Z'o|") - / , (1-23) 

n=l ' n=l 

where we have taken into account the completeness of the set |n). Obtained 
result is just the usual g.s. expectation value of an operator. 

The computation of static quantities (|1.22p still involves the summation 
over the complete set of Ngt states which is clearly not feasible in practice. 
To obtain a useful method, a further essential approximation replaces the 
full summation over \n) by a partial one over a much smaller set of random 
states [TTIHH]- Such an approximation is analogous to Monte Carlo methods 
and leads to a statistical error which can be well estimated and is generally 
quite small. Let us first consider only the expectation value (|1.19p with respect 
to a single random state |r), which is a linear combination of basis states 
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\r) = Y,Vrn\n), (1.24) 

n=l 

where rjm are assumed to be distributed randomly. Then the random quantity 
can be expressed as 

Ar = {r\e-^"A\r)/{r\e-^"\r) ^ 
= E V*rnVrm{n\e~^" A\m) / J2 €nr}rnMe'^"\m). (1.25) 

n,m=l ' n,m— 1 

Assuming that due to the random sign (phase) offdiagonal terms with ri^rjlrrm m ^ 
n on average cancel for large Ngt, we remain with 

^ = E \Vrnf {n\e~^"A\n) / ^ h™|'(n|e-'^^|n). (1.26) 

n—1 ' n—1 

We can express IryrnP = ^/Nst + Sm- Random deviations Sm should not be 
correlated with matrix elements 

therefore A^ is close to (A) with an statistical error related to the effective 
number of terms Z in the thermodynamic sum, i.e. 



Ar^ {A)il + Oil/VZ)), (1.27) 

Nst 

n n—1 

Note that for T oo we have Z — >■ Ngt and therefore at large Nst very 
accurate average (|1.28p can be obtained even from a single random state 
pTpS] . On the other hand, at finite T < oo the statistical error of Ar increases 
with decreasing Z. 

To reduce statistical error, in particular at modest T > 0, within the 
FTLM we sum in addition over R different randomly chosen \r), so that in 
the final application Eq. (|1.22l) leads to 

A/ 

r=l j=Q 

z^^Y.T.^~''^\m)\'- (1-29) 

Random states |r) = |(/)q) serve as initial functions for the Lanczos iteration, 
resulting in M eigenvalues with corresponding The relative statistical 
error is reduced by sampling (both for {A) and Z) and behaves as 

5{A)I{A) =0{II^/r1). (1.30) 
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For general operator A the calculation of \4>j) and corresponding matrix 
elements {tpj\A\r) is needed. On the other hand, the calculation effort is sig- 
nificantly reduced if A is conserved quantity, i.e., [-ff , ^] = 0, and can be 
diagonalized simultaneously with H. Then 

B. M 

- ^EE^"'''^'i(^i'/'Pi%- (1-31) 

r=l j"=0 

In this case the evaluation of eigenfunctions is not necessary since the element 
{r\ip^) = v^Q, Eq. (|1.6p . is obtained directly from eigenvectors of the tridiagonal 
matrix H^j. There are several quantities of interest which can be evaluated 
in this way, in particular the thermodynamic properties as internal energy, 
specific heat, entropy, as well as uniform susceptibility etc. [71ll9j. 

Taking into account all mentioned assumptions, the approximation (A) 
(|1.29p yields a good estimate of the thermodynamic average at all T. For low 
T the error is expected to be of the order of 0{1/\/R), while for high T the 
error is expected to scale even as 0{\/ \/NstR)- Since arguments leading to 
these estimates are not always easy to verify, it is essential to test the method 
for particular cases. 



1.4.2 Finite temperature Lanczos method: Dynamical response 

The essential advantage of the FTLM with respect to other methods is never- 
theless in the calculation of dynamical quantities. Let us consider the dynam- 
ical susceptibility as given by the autocorrelation function C(uj) (procedure 
for the general correlation function Cab{'^) is given in Ref. [7]), 

X"(t^) = 7r(l - e-P'^)C{uj), C{uj) = -Re / dte"^'C{t), (1.32) 

TT Jo 

with 

C{t) = (At(t)A(O)) = ^ Y,{n\e(-^+''^"A^e-'"'A\n). (1.33) 

n 

Expanding the exponentials in analogy to static quantities, EQ. (jl.20p . we get 
at) =Z-'J2f: (-^ + ^^)^ (-f)' (n|g^AtH'A|n). (1.34) 

n=l k.l=Q 

Expansion coefficients in Ea. (|1.34p can be again obtained via the Lanczos 
method, as discussed in Sec ll.4.11 Performing two Lanczos iterations with M 
steps, started with normalized |(/)q) — \n) and |^o) ^ ^l*^); respectively, we 
calculate coefficients Wki following the equation (|1.17p . We again note that 
(within the full basis |n)) the series are via Wki exactly evaluated within the 
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Lanczos basis up to order l,k < M. The latter yields through Eq. (|1.34p a 
combination of expansion, i.e. a combination of high-T, short-t (in fre- 
quency high-w) expansion to very high order. Extending and resuming series 
in k and I into exponentials, we get in analogy with Ea. (|1.20p 

N,t M 

at) = J2 E e-^''e''^''-'^Hn\^-)mA^i^-){i^-\A\n), (1.35) 

n— 1 

Finally replacing the full summation with the random sampling the FTLM 
recipe for the correlation function is 

R M 

C'M = E e-^^Hrm{rM^\r,){i'-\r}6{u;-i^^+el). (1.36) 

r— 1 i J — 1 

We check the nontrivial T — limit of above expression. If \n) are not 
orthogonal to the g.s. |!f'o)i then for large enough M the lowest-lying state 
converges to Eq ~ and iV'o) ~ I'^o), respectively. In this case we have 

R M 

C(c.,T = 0)«-^^^(^o|^^|^;)(i"|^|r)(r|tf'o)^(^ + i?o-£7) (1.37) 

At T '-^ one needs in general M ^ 100 in order that at least low-lying states 
relevant to Ea. (|1.37p approach |^") and — >■ Ej. Also a considerable 

sampling i? > 1 is required to get correct also amplitudes of separate peaks 
in the spectrum of Eq. (|1.37p which are a subject of a statistical error due to 
the incomplete projection on different random |r) in {ip"\A\r){r\WQ). Similar 
statistical error can in fact appear also for static quantities in Eq. (|1.29p . 



1.4.3 Finite temperature Lanczos method: Implementation 

Most straightforward is the implementation of the FTLM for static quantities, 
Eq. (|1.29p . In particular for conserved quantities, Eq. ()1.3ip . the computation 
load is essentially that of the g.s. Lanczos iteration repeated R times and only 
a minor changes are needed within the usual g.s. Lanczos code. 

On the other hand, for the dynamical correlation function (|1.36p the 
memory requirement as well as the CPU time is dominated mostly by the 
evaluation of the matrix element (V'fl^^lV'J) where the operations scale as 
Op cx RM'^Nst and memory as Mem oc MNgf This effectively limits the ap- 
plication of the FTLM to 50 < M < 500 where the lower bound is determined 
by the convergence of the g.s. jtf'o)- Still, it should be noted that the calcula- 
tion can be done simultaneously (without any additional cost) for all desired 
T, since matrix elements are evaluated only once. Evidently, one should use 
as much as possible symmetries of the Hamiltonian, e.g., A^e, •S'foj, q to reduce 
effective N^t by splitting the sampling over different symmetry sectors. 
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The effect of finite M is fess evident. Since M ^ 100 is enough to converge 
well few lowest levels, it is also generally satisfactory for reliable dynamical 
correlation functions at low T. At high T, however, one can observe very regu- 
lar oscillations which are the artifact of the Lanczos iterations with M <C Ngt ■ 
Namely, the procedure generates between extreme eigenvalues quite equidis- 
tant spectrum of quasi-states with the level spacing Ae ~ AE/M, where AE 
is the full energy span of MB eigenstates. The effect is well visible in Fig ll.2l 
where the high-T result for the spin structure factor S{q = 7t,uj) for the ID 
Heisenberg model, Ea. (|l.ip . is presented for various M. It is evident that for 
the presented case (TV — 24 and AE ~ 16 J) M > 200 is sufficient to obtain 
smooth spectra even for high T ^ J. However, larger M are advisable if 
sharper structures persist at high T. 




Fig. 1.2. High-T spin structure factor S{q — 7r,tj) for the ID Heisenberg model, as 
calculated with different number of Lanczos steps M. 



The role of random sampling R is less important for intermediate and 
high T since the relative error is largely determined via Z as evident from 
Eq. (|1.28p . Larger i? ^ 1 is necessary only for the correct limit T — >■ (for 
given system size) and for off-diagonal operators A. 

One can claim that the FTLM in general obtains for all reachable systems 
results which are at any T very close to exact (full ED) results for the same 
finite (given N) system and the accuracy can be improved by increasing M 
and R. Still, it remains nontrivial but crucial to understand and have in control 
finite size effects. 

At r = both static and dynamical quantities are calculated from the 
g.s. |!fb)i which can be quite dependent on the size and on the shape of the 
system. At least in ID for static quantities the finite-size scaling N ^ oo can 
be performed in a controlled way, although in this case more powerful meth- 
ods as, e.g., the DMRG are mostly available. In higher dimensional lattices. 
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e.g., in 2D systems finite-size scaling is practically impossible due to very re- 
stricted choice of small sizes and different shapes. Also g.s. {T — 0) dynamical 
quantities are often dominated by few (typically Np < M) peaks which are 
finite-size dominated On the other hand, T > generally introduces the 
thermodynamic averaging over a large number of eigenstates. This directly 
reduces finite-size effects for static quantities, whereas for dynamical quanti- 
ties spectra become denser. From Eq. (jl.36|) it follows that we get in spectra 
at elevated T > typically Np cx RM^ different peaks resulting in nearly 
continuous spectra. This is also evident from the example of a high-T result 
in FigO 

It is plausible that finite-size effects at T > become weaker. However, it 
should be recognized that there could exist several characteristic length scales 
in the considered physical (and model) system, e.g. the antiferromagnetic 
(AFM) correlation length ^, the transport mean free path Ig etc. These lengths 
generally decrease with increasing T and results for related quantities get a 
macroscopic relevance provided that $,{T),ls{T) < L where L cx N^^'^ is the 
linear size of the system. However, there exist also anomalous cases, e.g., in 
an integrable system Is can remain infinite even at T — )■ oo [201121] . 
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Fig. 1.3. Finite-size temperature Tfs vs. hole doping Ch in the 2D t-J model with 
J/t = 0.3, as calculated with the FTLM in system of iV = 18 sites [7]. 

A simple criterion for finite size effects one can use the normalized ther- 
modynamic sum Z{T), Ea. (|1.28p . which provides the effective number of MB 
states contributing at chosen T (note that for a system with a nondegenerate 
g.s. Z{T = 0) = 1). Finite-size temperature Tfs can be thus defined with the 
relation Z{Tfs) = Z* where in practice the range 10 < Z* < 50 is reasonable. 
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Clearly, the FTLM is best suited just for systems with a large density of low 
lying MB states, i.e., for large Z at low T . 

Since ^(T) is directly related to the entropy density s and the specific heat 
Ct, of the system, large Z at low T is the signature of frustrated quantum MB 
systems which are generally difficult to cope with other methods (e.g., the 
QMC method). Such are typically examples of strongly correlated electrons 
with an inherent frustration, e.g., the doped AFM and the t-J model, Eq. (|1.2p . 
in the strong correlation regime J < t. We present in Fig ll.3l as an example 
the variation of Tfs in the 2D t-J model with the hole doping Ch = Nh/N, as 
calculated for different Z* — 30 — 300 for the fixed system of iV = 18 sites and 
J/t — 0.3 as relevant for high-T cuprates. It is indicative that Tfs reaches the 
minimum for intermediate (optimum) doping = cj^ ~ 0.15, where we are 
able to reach Tfg/t ~ 0.1. Away from the optimum doping Tfg is larger, i.e., 
low-energy spectra are quite sparse both for undoped AFM and even more 
for effectively noninteracting electrons far away from half-filling (for nearly 
empty or full band). 

1.4.4 Low Temperature Lanczos Method 

The standard FTLM suffers at T — from a statistical error due to fi- 
nite sampling R, both for the static quantities, Eqs. (|1.29|) . (11.301) . as well 
as for the dynamical correlations, Eas. (|1.36p . (|1.37p . The discrepancy can be 
easily monitored by the direct comparison with the g.s. Lanczos method, 
Eas. (|1.8ll.lip . To avoid this problem, a variation of the FTLM method, named 
Low-temperature Lanczos method (LTLM) has been proposed [5] which ob- 
tains correct g.s. result (for finite systems) independent of the sampling R. 
The idea of LTLM is to rewrite Eq. (|1.19p in a symmetric form 

{A) = ^j:{n\e-^"/'Ae-^"/M, (1.38) 

n=l 

and insert the Lanczos basis in analogy with the FTLM, Eq. ()1.19| ). now rep- 
resented with a double sum 

i? A/ 

= E e-^^^'^+<'>^'{r\r,){r,\m){^l\r), (1.39) 

The advantage of the latter form is that it satisfies the correct T = limit 
provided that the g.s. is well converged, i.e., |i/'o) ^ \^o)- It then follows from 
Eg. lfOa . 

R. R 

{A) = Y,{rWo){MAm{Mr)/Y.^Mr){rWo) = ('^ol^l'^o), (1-40) 

r— 1 r— 1 

for any chosen set of |r). For the dynamical correlations C{t) one can in 
straightforward way derive the corresponding expression in the Lanczos basis 
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r=l 

X (5(.;-6-' + i(6[ + 6[)). (1.41) 

It is again evident that for T — >■ the sampling does not influence results being 
correct even for i? = 1 if the g.s. \'Fq) is well converged for all starting \r). 
The payoff is in an additional summation over the new Lanczos basis which 
starts from each AjV'f) in Eq. (|1.40l) . Since the LTLM is designed for lower 
T there one can effectively restrict summations in (i,/) in Ea. (jl.4ip to much 
smaller M' <C M where only lowest states with e[ , e[ ~ Eq contribute [5] , and 
in addition use smaller Mi <C M for the basis • 

An alternative version for Lanczos-type approach |22] for dynamical quan- 
tities is not to start the second Lanczos run from A\r) [7] or from [5], 
but from 

M 

\Ar}=J2m)e-^''^'{m. (1.42) 

1=0 

In this way one obtains with the second Lanczos run the Lanczos eigenstates 
\il'k)j which cover the relevant Hilbert space for starting random vector |r) and 
the inverse temperature /3. The resulting dynamical autocorrelation function 
is 

r— 1 i,k—0 

X Siuj - + £[)• (1-43) 

In this way the sufficiency of only one random vector in the T — limit is 
reproduced, while at T > the algorithm has the same time efficiency as the 
FTLM, but with much smaller random sampling needed to reach the same 
accuracy (at least for low T). However, the price paid is that results for each 
T need to be calculated separately, while within the FTLM all T (or T up to 
certain value within the LTLM) are evaluated simultaneously. 



1.4.5 Microcanonical Lanczos Method 



While most investigations in strongly correlated systems focus on the low- 
T regime, there are systems where dynamical properties are nontrivial even 
at high T. Well known such case is the spin diffusion constant Ds{T) in 
the isotropic Heisenberg model, Eq. ()l.l|) . which is not known by value and 
moreover not even its existence at any T > 0. Similar although somewhat less 
controversial is the case of transport quantities, both for integrable or generic 
nonintegrable models. Whereas the FTLM seems well adapted for studies 
of transport response functions, oscillations due to limited M can affect the 
crucial low-w resolution as seen also in Fig ll.2l 
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At elevated T it is therefore an advantage to use the Microcanonical Lanc- 
zos method (MCLM) [5], employing the fact from statistical physics that 
in the thermodynamic limit (for large system) the microcanonical ensemble 
should yield the same results as the canonical one. The shortcoming of the 
MCLM emerges since in finite systems statistical fluctuations are much larger 
within the microcanonical ensemble. Still, reachable finite-size systems have 
very high density of states in the core of the MB spectrum as probed by high 
T. Hence, statistical fluctuations are at high T effectively smoothed out in 
contrast to low-T properties dominated by a small number of low lying MB 
states. 

The implementation of the MCLM is quite simple and straightforward. 
One first determines the target energy A = {H)(T) which represents the mi- 
crocanonical energy equivalent to the canonical average energy for chosen T 
and the system size N. Since A is a parameter within the MCLM, one can re- 
late it to T by performing either FTLM (simplified due to conserved quantity 
H) on the same system or extrapolating full ED results (with linear depen- 
dence on N) on small lattices. Next we find a representative microcanonical 
state for the energy A. One convenient way within the Lanczos-type ap- 
proach is to use the new operator 

V = {H-Xf. (1.44) 

Performing Lanczos iterations with the operator V yields again the extremum 
eigenvalues, in particular the lowest one close to V ^ 0. In contrast to the 
g.s. procedure, the convergence to a true eigenstate cannot be reached in 
system sizes of interest even with Mi ^ 100. The reason is extremely small 
eigenvalue spacing of operator V scaling as AVn oc {AE/Ngt)"^, AE being the 
whole energy span within the given system. Fortunately such a convergence is 
not necessary (even not desired) since the essential parameter is small energy 
uncertainty ge, given by 

al = {^x\V\^x). (1.45) 

For small energy spread ge/AE < 10~^ typically Mi ~ 1000 is needed. 
Again, to avoid storing Mi Lanczos wavefunctions Lanczos procedure is 
performed twice as described in Sec ll.2.2] i.e., the second time with known 
tridiagonal matrix elements to calculate finally I'l^x) in analogy with Eq. (|1.6p . 
The latter is then used to evaluate any static expectation average (A) or the 
dynamical correlation function as in Ea. (ll.9p . 

C{io, A) = {^x\A^ \ mx). (1.46) 

The latter is evaluated again using Lanczos iterations with M2 steps starting 
with the initial wavefunction |0o) oc and C{uj) is reprepresented in terms 

of continued fractions. Since the MB levels are very dense and correlation 
functions smooth at T ^ large M2 ^ 100 are needed but as well easily 
reachable to achieve high-w resolution in C{uj). 
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It is evident that the computer requirement for the MCLM both regarding 
the CPU and memory are essentially the same as for the g.s. dynamical cal- 
culations except that typically Mi, M2 ^ 100. In particular requirements are 
less demanding than using the FTLM with M > 100. A general experience is 
that for systems with large Nst >• 10000 the MCLM dynamical results agree 
very well with FTLM results for the same system. It should be also noted 
that actual frequency resolution Sui in C{uj), Eq. (|1.46|) . is limited by Stu ^ ue 
which is, however, straightforward to improve by increasing Mi, M2 with typ- 
ical values Ml, M2 > 1000. One can also improve MCLM results for any T by 
performing an additional sampling over initial random starting |(/)o) as well as 
over A with a probability distribution p{X) simulating the canonical ensemble 
in a finite-size system, i.e., by replacing Ea. (|1.46p with 



1.4.6 Statical and dynamical quantities at T > 0: Applications 

The FTLM has been designed to deal with simplest tight-binding models of 
strongly correlated electrons, at the time mostly with challenging microscopic 
electronic models of high- Tc superconductors [6l[7], where besides supercon- 
ductivity there is a variety of anomalous non- Fermi-like properties even in the 
normal state. Clearly of interest in this connection are prototype MB models 
as the Heisenberg model, Eq. (|l.ll) . the t-J model, Eq. (|1.2l) and the Hubbard 
model on the 2D square lattice. Unfrustrated Heisenberg model can be numer- 
ically studied on much bigger lattices with QMC and related methods. The 2D 
Hubbard model was and still is mostly subject of DMFT and QMC studies, 
since at half-filling or close to it the Lanczos methods are quite restricted due 
to large Nat even for modest sizes N ^ 16. Therefore one focus of Lanczos- 
based approaches was on the t-J model being with some generalizations a 
microscopic representation of electronic properties of high- Tc cuprates. 

Thermodynamic quantities as chemical potential /x, entropy density s, spe- 
cific heat Cy are the easiest to implement within the FTLM. Their T- and 
(hole) doping c?i-dependence within the t-J model on a 2D square lattice (cal- 
culated up to TV = 26 sites) reveal already very anomalous behavior of doped 
Mott insulators [23] (as evident already from Fig |1.3[) . confirmed also by re- 
sults for more complete Hubbard model [24] . An introduction of next-neighbor 
hopping t' generates also an asymmetry in thermodynamic properties between 
hole-doped and electron-doped cuprates [33] as well quite dramatic influence 
of stripe order [55] consistent with the physics of cuprates. 

The advantages of the FTLM and also its feasibility for the 2D t-J model 
are even more evident in quite numerous studies of spin and charge dynamics 
at T > [7] which show good agreement with neutron scattering and NMR 
|19I27IB5] . optical conductivity cr(a;) and resistivity p{T) [30,31], Hall constant 
Rh{T) [32] and a general non- Fermi- liquid behavior of cuprates [5^, as well as 




(1.47) 
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the puzzling strong influence of nonmagnetic impurities [33^ . As an example 
of a transport quantity hardly accessible by other methods we present in 
Fig ll.4l the universal planar resistivity p{T), as extracted from the dynamical 
conductivity 17(0; — >■ 0) = 1/p, within the t-J model for different doping 
levels Ch [31]. Result in Fig ll.4l clearlv shows a linear dependence below the 
pseudogap temperature T* dependent on doping Ch- Another characteristic 
signature is a saturation (plateau) of p{T) at low doping and the universal 
trend at high T. 
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Fig. 1.4. Normalized 2D resistivity Chp vs. T/t within the t-J model with J /t = 0.3 
for different hole concentrations Ch 1311 . 



Spectral properties as manifested in single-particle spectral functions 
^(k, w) are in the core of the understanding of cuprates, as well as of strongly 
correlated electrons in general. Here, even g.s. and low-T properties are the 
challenge for numerical studies whereby the FTLM can be viewed as a con- 
trolled way to get reliable (macroscopic- like) T — ^ result, in contrast to 
quite finite-size plagued results obtained via g.s. Lanczos procedure [5j. Using 
the FTLM at T ~ Tfs with the twisted boundary condition can simulate a 
continuous wavevector k. Using in addition the course graining averaging one 
can reach results for A(k, oj) giving insight into the electron vs. hole 

doped angle-resolved photoemission experiments, quasiparticle relaxation and 
waterfall-like effects. A characteristic result of such studies is in Fig ll.Sl for the 
single-particle density of states N{ijj) = A(k, w) [34|. Here, the strength 
of the FTLM is visible in the high uj resolution within the most interesting 
low-w window. Interesting and reproducible are also nontrivial spectral shapes 
as the sharp peak close to oj < and a broad shoulder for cj ^ 0. Most im- 
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portant is, however, the evident pseudogap (observed also experimentaUy in 
cuprates) visible at a; ^ in the low-doping regime. 




Fig. 1.5. Density of states M{ui) for different dopings Ch within the extended t-J 
model with n.n.n. hopping t' — ~0.3t and t' = 0, respectively [34) . 

Besides the challenging models for cuprates there have been also studies of 
static and dynamical properties of multiband and multiorbital models which 
either reduce to the generalized t-J model [37] or to Kondo lattice models 
[55Wn] and the Falicov-Kimball model [H]. While the increasing number of 
local basis states K clearly limits the applicability of ED-based methods, they 
are competitive in treating nontrivial frustrated spin models less suitable for 
the QMC and other methods, however closely related to physics of novel 
materials. Moreover, frustrated models are characterized by a large entropy 
density s and related low Tfs essential conditions for feasibility of FTLM 
results. Examples of such systems are the Shastry-Sutherland model [42ll43) . 
2D J1-J2 model [44], and properties of frustrated magnetic molecules [45U47) . 

Another class of problems which can be quite effectively dealt with the 
FTLM and MCLM approaches is the fundamental as well as experimentally 
relevant problem of transport in ID systems of interacting fermions as realized, 
e.g., in quasi-lD spin-chain materials [48]. It has been recognized that the 
transport response at any T > crucially differs between integrable and 
nonintegrable systems. Since the ID isotropic as well as anisotropic Heisenberg 
model, Eg. dl.ip . is integrable it opens a variety of fundamental questions 
of anomalous transport in such systems, the effects of perturbative terms 
and impurities. Such fundamental questions on transport and low-w dynamic 
response remain nontrivial even at high T [20l[2T], hence the MCLM is the 
most feasible and straightforward method. It has been in fact first probed on 
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the anomalous transport in ID insulators [35^ but furtheron used to study 
interaction- induced transport at T > in disordered ID systems [50ll51) . 
incoherent transport induced by a single either static [5 2) or dynamical spin 
impurity [53 . 

In Fig ll.61 we present as an example the MCLM result for the dynam- 
ical spin conductivity in the anisotropic Heisenberg model, Ea. (|l.ip . where 
J^^ 7^ J^^ = = J in the Ising-like (with the spin gap in the g.s.) 
regime A = J^^/J > 1. Results for the high-T dynamical spin conductivity 
Tij{lo) are shown for various next-neighbor (anisotropic) coupling a = J|^/ J. 
First message is that the MCLM as the method is well adapted for the high- 
w resolution (here using Mi = = 2000) and reaching large TV = 30 
{Nst ~ 5.10^ in a single — 0,q sector). Another conclusion is that the 
dynamics of such systems is very anomalous. For the integrable case a = we 
find (To = cr(0) ^ but also an anomalous finite-size peak at ujp cx 1/7V 09]. At 
the same time breaking integrability with a > appears to lead to <to > still 
approaching an 'ideal' insulator (insulating at all T) for a weak perturbation 
(7o(a-^ 0) ^ [M]. 



1.5 Reduced Basis Lanczos Methods 

The main shortcoming of ED approaches are finite-size effects that tamper 
calculations on small lattice systems. Exponentially growing Hilbert spaces 
represent the main obstacle against extending ED calculations to larger lat- 
tices. One way to extend ED calculations is to reduce the complete Hilbert 
space and keep only states that give a significant weight in the g.s. wavefunc- 
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tion or in the relevant Hilbert space, e.g., at T > 0. Here, the crucial step 
represents developing a feasible algorithm for the basis reduction. 

One clear example of very effective basis construction is the DMRG 
method, most feasible for ID correlated systems. Within this method, de- 
scribed in other chapters, some intermediate steps of ED diagonalization are 
frequently also performed via Lanczos iterations [55'. Moreover, there are 
recently developments which are extending the application of DMRG proce- 
dure to T > dynamical response [56l[57] involving more directly features 
of Lanczos-based approaches. In particular, recent finite-T dynamical DMRG 
(FTD-DMRG) method [S^ combines the DMRG selection of MB states and 
the LTLM procedure as described in Sec fTXSl and in Eqs. ([L42| . ([L43)) to 
evaluate dynamical correlations at T > 0, both being effective at low T. So 
far the FTD-DMRG method has been applied to find novel features within 
the ID J1-J2 model [SS]. 

The ideas to find the reduced basis sets is, however, more general going 
beyond the ID systems and has been successful in solving several problems 
of strongly correlated systems. Even before the discovery of high- Tc super- 
conductors that boosted research on models of correlated systems, Brinkman 
and Rice |59j developed a string representation of configurational space to 
compute the band renormalization and mobility in the atomic limit of the 
Hubbard model. The string picture was later used to compute single hole 
properties and to estimate hole-pair binding within the t-J model by many 
authors [60] . In connection with ED calculations on finite lattices the trunca- 
tion scheme based on the limitation of the maximal length of strings leads to 
slow convergence in terms of the number of reduced basis states |61| in the 
spin-isotropic limit of the t-J model. Rather slow convergence is achieved also 
using a cluster diagonalization scheme together with the systematic expansion 
of the Hilbert space [62] . 

The exact diagonalization method in limited functional space (EDLFS) 
was originally developed specifically for solving the problem of a single charge 
carrier (hole) doped in the AFM background [S3]. The basic principle of the 
method is based on the above mentioned string picture, that emerges as a 
moving hole through the AFM background creates in its wake paths of over- 
turned spins - strings. EDLFS was later generalized to include phonon degrees 
of freedom [64][65] . For the purpose of presenting the method we start by con- 
sidering the prototype t-J model, Ea. (|1.2p . coupled to dispersionless optical 
phonons on a square lattice 

H = Ht+ Hj + Hg + = 

+ 5E(^""*)("i^ + ''*)+'^oE^^"*' (^-"^^^ 

i i 
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where t represents the nearest-neighbor hopping, is the phonon annihilation 
operator, Ui = riis and all double sums run over the nearest neighbor pairs. 
The third term represents electron-phonon coupling g and the last the energy 
of Eistein-type phonons ujq . 

• ••••• 




• ••••• 

A B 

Fig. 1.7. (Color online) Schematic representation of a particular basis wavefunction 
obtained using — 10, m — 3, and A^;, — 3. Circle represents the hole position, 
crosses portray spin flips and numbers indicated excited phonon quanta. Dots repre- 
sent lattice sites with no spin flips and zero excited phonon quanta. In this particular 
case, Nfi = 6 and Np = 4. Presented basis function is one of a total of N^t = 42 x 10^ 
states, generated using Ea. (|f .49p . 



The Neel state is the g.s. of the model when t = j — g = 0. We construct 
the limited functional space (LPS) for one hole starting from a Neel state on 
an infinite 2D lattice with a hole placed in the origin and zero phonon degrees 
of freedom \(f)o) = co|Neel;0). We then proceed with the generation of new 
states, using the following generator 

{m = {m+Y,Hir'\M, (1.49) 

i=l 

where Ht and Hg are off-diagonal terms in Eq ll.48l respectively. We se- 
lect Us — l,Ns and a fixed value of m. This procedure generates strings 
with maximum lengths given by Ns and phonon quanta at a maximal dis- 
tance A'^s — 1 from the hole. Parameter m provides the creation of additional 
phonon quanta leading to a maximum number of phonons Nph = mNg in 
the system. Typically, larger Np^ are necessary to achieve convergence in the 
strong electron-phonon coupling (polaron) regime. While constructing LPS 
we take into account explicitly translational symmetry. Le., we store only one 
out of (in principle) infinitely many translationally invariant states, called a 
parent state. However, the Neel state is a state with a broken translational 
symmetry so allowed translations are generated by two minimal translations 
ri^2 = (1,±1). The basis wavefunctions are represented by one of the two 
nonequivalent hole positions r^, sets of strings, representing overturned spins 
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relative to the Neel state (spin flips), and occupation numbers representing 
excited plionon quanta, 

\4>) = |r/j;ri,r2,...,rAr^.,;nr;,"-r;,,---,"-r'„J, (1-50) 

where represent spin-flip coordinates, n-r' number of phonon quanta at site 
and Nfi < Ng is the total number of spin flips. 

Since the Hilbert space grows exponentially with Ng there exist many 
intuitive physically motivated restrictions optimizing LFS's for different pa- 
rameter regimes to slow down the exponential growth. A systematic control 
over the growing Hilbert space can be achieved through realization that the 
diagonal energy of long strings grows linearly with Ns. Assuming that the 
contribution of long strings in the g.s. wave-function is negligible, we in- 
troduce an additional parameter Nb < Ns that restricts generation of long 
strings by imposing a condition under which all coordinates of spin flips sat- 
isfy l/x/j — < Nb;fi — {x,y} where h and / refer to hole and spin-flip 
indexes, respectively. Application of this condition improves the quality of 
the LFS by increasing the number of states containing spin flips in the vicin- 
ity of the hole while keeping the total amount of states within computationally 
accessible limits. 

Fig ll. 71 represents a particular state generated using Ea. (|1.49l) . Black and 
grey dots represent sites with spins 'up' and 'down', respectively. The hole 
starts at the position (—1,1) in the direction (0,-1), and travels along the 
path indicated by arrows. Crosses represent spin-flips and numbers represent 
excited phonon quanta, generated along the hole's path. The effect of the 
parameter m is to set the maximal number of excited phonon quanta on 
indicated positions. In this particular parent the hole ends its path on the 
B-sublattice. Using allowed minimal translations, the position of the hole in 
the parent state is at Vh = (1, 0). 

After the generation of LFS the full Hamiltonian in Eq. (|1.48|) is diagonal- 
ized within this LFS using the standard Lanczos procedure. The efficiency of 
the EDLFS method in case of the J model and stability of results against 
varying Ng and Nb has been shown in detail in Ref. 63j . 

Besides reliable results obtained using the EDLFS there are other im- 
portant advantages over most other methods: a) it is highly efficient, b) the 
method is free of flnite-size effects, c) it treats spin and lattice degrees of 
freedom on the same footing while preserving the full quantum nature of the 
problem, and d) it allows for computation of physical properties at an arbi- 
trary wavevector. Even though results depend on the choice of parameters 
deflning the functional generator in Eq. (|1.49|) . such as Ng and m, reliable re- 
sults are obtained already for relatively small sizes of the Hilbert space Ngt , 
typically up to three orders of magnitude smaller than in the case of ex- 
act diagonalization techniques. For most static as well as dynamic quantities 
convergence to the thermodynamic limit can be achieved with a systematic 
increase of Ng and m. 
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The EDLFS obviously has few Hmitations: a) the method is hmited to 
calculations in the zero-doping limit, e.g., = 1,2 mobile particles immersed 
in an infinite (AFM) background, b) the spin anisotropy is inherently built in 
the method, and c) due to the broken translational symmetry of the starting 
wavefunction, calculations are limited to the reduced AFM Brillouin zone. 

Due to its high efHciency in dealing with spin fluctuation the method 
represents one of the few successful methods that allows the addition of lattice 
degrees of freedom to a correlated electron model. The EDLFS handles spin 
and lattice degrees of freedom by preserving the full quantum mechanical 
nature of the problem and enables a direct calculation of the dynamic response 
functions in terms of real frequency uj. 

The EDLFS can be rather in a straightforward way generalized for the 
study of the two- hole = 2 problem ^66,t67j. ED studies in the t-J model 
for Nh — 2 were performed on a 2D square lattices up to = 32 sites [55]. 
Still, the maximal distance between two holes remains rather small Z,nax = 
yj N/2 = 4. Attempts to increase the lattices sizes beyond ED studies led 
authors to investigate various truncated basis approaches [6T1|62] and using 
small sets of variational wavefunctions with given rotational symmetries based 
on 'string' and 'spin - bag' pictures [69] . 




While the problem of two holes within the t-J model represents a challeng- 
ing problem, the addition of quantum phonons seems an almost unachievable 
task. Nevertheless, the EDLFS due to an efficient choice of LFS can handle 
this problem well. The construction of the LFS starts from a Neel state with 
holes located on neighboring sites and with zero phonon quanta. Such a state 
represents a parent state of a translationally invariant state. We generate new 
parent states in analogy with Eq. (|1.49|) by applying the generator of states 



where Hj denotes the off-diagonal spin exchange term in Ea. (jl.48|) applied 
only to erase spin flips generated through application of Ht. This allows the 
creation of states with holes positioned further apart that are not connected 
with spin strings. Note that for larger Ns > Q some of such states would be 
generated even without evoking Hj term via the Trugman loops [BD] • One of 
the advantages of this method in comparison to other approaches is that it 
allows much larger distances between the holes, /max = A^s + 1- Note also that 
the functional generator preserves as well the point-group symmetry. 

As an example we present in Fig ll.Sl the probability of finding a hole pair 
at a distance r 



where nj* = 1 — is local hole density. In the case of dimensionless electron- 
phonon coupling A — g'^/8u!ot — 0, the maximal allowed distance between 





(1.51) 




(1.52) 
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Fig. 1.8. Probability P{r) of finding a hole pair at a distance r within the t- 
J-Holstein model with J/t — 0.3 and iOo/t = 0.2, evaluated for electron-phonon 
coupling parameters A = and A — 0.25. Different Ns and m as indicated in legends 
were chosen giving Nst ~ 7.f0'' and Nst ^ 26. fO®, respectively. Both solutions 
correspond to a d-wave symmetry. 



holes is Zmax = 13 while for A = 0.25 Zmax = 7, whereby results for A = are 
consistent with previous findings using ED [68] , reduced-basis ED [62] , as well 
as string picture variational approaches [69] . 



1.6 Real Time Dynamics using Lanczos Method 

Research in the field of non-equilibrium dynamics of complex quantum sys- 
tems constitutes a formidable theoretical challenge. When dealing with ED 
approaches or calculations in the reduced basis, the time evolution of the time 
- dependent Shrodinger equation, 

^mp_^^^^^^^^^^ (1.53) 

can be efficiently obtained using the time - dependent Lanczos technique, 
as originally described in Ref. |70| and later applied and analysed in more 
detail [71]. One of the straightforward reasons is that most commonly the 
Lanczos method is used to compute g.s. of MB Hamiltonian. Generalizing 
the method to time - dependent calculation represents only a minor change 
to already existing codes. Even though the method is most suitable for the 
time evolution of the time - independent Hamiltonian, it can nevertheless be 
applied even to the time - dependent case. The time evolution of is 
then calculated by step-vise change of time t in small time increments St, 
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generating at each step Lanczos basis of dimension M (typically M < 10), to 
follow the evolution 

M 

\^{t + St)) e-*^«**|*'(0) ^ J2 e-"''*IV'i)(V';|tf'(i)), (1-54) 

1=1 

where \'ipi),ei,l — 0, M are Lanczos eigenfunctions and eigenvalues, respec- 
tively, obtained via the Lanczos iteration started with |0o) — l'^(^))- The ad- 
vantage of the time-evolution method following Eq. (|1.54p is that it preserves 
the normalization of \W{t + St) for arbitrary large dt. The approximation of fi- 
nite M in Ea. (|1.54p is also correct at least to the M-th Taylor-expansion order 
in St. It is, however, important to stress that St should be chosen small enough 
to take into account properly the time-dependence of H{t). E.g., in the case 
of driving the system with an external constant electric field, St/ts ^ 10~^ 
where is the Bloch oscillation period [54 1I72H74] . 

So far, investigations of correlated driven systems under the influence of a 
driving electric field in ID systems using the Lanczos time evolution method 
focused on generic systems, like the metallic and Mott - insulating regime 
within the ID model of interacting spinless fermions [54l[72]. Even though 
rather small systems can be studied it has been established that steady state 
can be reached without any additional coupling to a heat bath, provided that 
the Joule heating of the system is properly taken into account. 

The case of a single charge carrier in an inelastic medium driven by 
the external electric field in ID as well as 2D has been investigated using 
EDLFS combined with a Lanczos-based time-evolution approach [72 l [74 j [75 ] . 
The strength of EDLFS is in construction of the Hilbert space that enables 
not only an accurate description of the ground state of the single carrier sys- 
tem but it allows for enough extra inelastic (spin or phonon) excitations to 
absorb energy, emitted by the field driven carrier, until the system reaches 
the steady state. This again enables a proper description of the steady state 
without coupling the system to an external thermal bath. 

1.7 Discussion 

Exact diagonalization based methods, both the full ED and Lanczos-type ED 
approach, are very extensively employed in the investigations of strongly cor- 
related MB quantum systems in solid state and elsewhere. The reason for their 
widespread use are several: a) unbiased approach to the MB problem without 
any simplifications or approximations, independent of complexity of the MB 
system, b) relative simplicity of generating the codes for various models and 
observables, c) easy and straightforward testing of codes, d) direct interpre- 
tation of obtained quantum MB states and their possible anomalous struc- 
ture and properties, e) high pedagogical impact as a quick and at the same 
time very nontrivial introduction into the nature of MB quantum physics. 
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Also the Lanczos-based methods described in this review, i.e., the g.s. Lanc- 
zos method for static and dynamic quantities, and somcwliat more elaborate 
FTLM, MCLM, LTLM and EDLFS, require rather modest programming ef- 
forts in comparison with more complex numerical methods, e.g., the QMC- 
and DMRG- based methods, as described in other chapters. 

Clearly, the main drawback of ED methods is the smallness of lattice sizes 
N determined by a limited number of basis states (at present Nst < 10^) 
treated with a Lanczos iteration procedure. The achievable N with ED meth- 
ods appears quite modest in comparison with some established and recently 
developed numerical methods, as the QMC, DMRG, matrix-product-states 
methods etc. Still, in spite of very intensive development and advance of novel 
numerical methods in last two decades there are several aspects of strong- 
correlation physics, where ED-based methods are so far either the only feasible 
or at least superior to other methods. In this chapter we have focused mostly 
on Lanczos-based methods and applications where they are still competitive 
and get nontrivial results with a macroscopic validity: 

a) MB g.s. and its properties of frustrated and complex models mostly so far do 
not offer alternative powerful methods except ED approaches, at least beyond 
D = 1 systems where DMRG-based methods can be effectively applied, 

b) T > static properties evaluated with Lanczos-based methods as the 
FTLM and the LTLM are as well most powerful and reliable for frustrated 
and complex system, in particular in systems with high degeneracies of MB 
states and large entropy at low T, 

c) T > Lanczos methods for dynamical quantities, as the FTLM and MCLM, 
yield for many models and geometries results superior to other methods or 
even the only accessible results in several cases. In particular the advantage of 
the latter methods is high ui resolution at all temperatures beyond the finite 
size limit T > Tfs, the macroscopic- like results at low T with a proper scaling 
to T ^ 0, and the possibility of detailed studies of systems with nontrivial 
(anomalous) dynamics at any, in particular high T. 

d) Lanczos technique of ED is also the natural application within the methods 

with a restricted MB basis sets as the EDLFS and DMRG-typc; targeting as 
well as in the real-time evolution studies of strongly correlated systems. 
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